function model_verify(model_settings,data,options)
h = zeros(length(data(1).t),1);
rho    = h;
meanR2 = h;
varR2  = h;
for i = 1:length(model_settings.W)
    t = data(i).t;
    dt = t(2)-t(1);
    h(1) = 0;
    rho(1) = 0;
    meanR2(1) = 0;
    varR2(1)  = 0;
    State.h   = 0;
    State.rho = 0;
    State.meanAlpha2 = zeros(1,20);
    State.meanBeta2  = zeros(1,20);
    State.meanR2     = 0;
    State.varAlpha2  = zeros(1,20);
    State.varBeta2   = zeros(1,20);
    State.varR2      = 0;
    for j = 2:length(t)
        State   = model(model_settings.W(i),State,dt,model_settings);
        h(j)      = State.h;
        rho(j)    = State.rho;
        meanR2(j) = State.meanR2;
        varR2(j)  = State.varR2;
    end
    msg_title = sprintf('w = %f',model_settings.W(i));
    if ~isempty(data(i).h)
        figure;
        subplot(2,1,1);
        plot(t,data(i).h,t,h);
        yalbel('h');
        xlim([0,options.t]);
        title(msg_title);
        subplot(2,1,2);
        plot(t,data(i).rho,t,rho);
        ylabel('SOR');
        xlim([0,options.t]);
    end
    figure;
    subplot(2,1,1);
    plot(t,data(i).meanR2,t,meanR2);
    xlim([0,options.t]);
    ylabel('mean(R2)');
    title(msg_title);
    subplot(2,1,2);
    plot(t,data(i).varR2,t,varR2);
    ylabel('var(R2)');   
    xlim([0,options.t]);
end